Quantum Entangled Dark Solitons Formed by Ultracold Atoms in Optical Lattices 
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Inspired by experiments on Bose-Einstein condensates in optical lattices, we study the quantum 
evolution of dark soliton initial conditions in the context of the Bose-Hubbard Hamiltonian. An 
extensive set of quantum measures is utilized in our analysis, including von Neumann and generalized 
quantum entropies, quantum depletion, and the pair correlation function. We find that quantum 
effects cause the soliton to fill in. Moreover, soliton-soliton collisions become inelastic, in strong 
contrast to the predictions of mean-field theory. These features show that the lifetime and collision 
properties of dark solitons in optical lattices provide clear signals of quantum effects. 
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Systems of ultracold atoms loaded into optical lat- 
tices offer an excellent experimental and theoretical test 
bed for the study of complex quantum many-body phe- 
nomena, including out-of-cquilibrium quantum dynam- 
ics. The unprecedented tunability of system parame- 
ters such as lattice height, filling, dimensionality, and 
species allows one to experimentally simulate key con- 
densed matter Hamiltonians This control also per- 
mits study of dynamical properties of the system, an 
aspect of quantum lattice physics not typically accessi- 
ble in the solid state. Examples include the experiment 
of the nonadiabatic transition across the Mott-superfluid 
border (2] , the experiment of strongly damped quantum 
transport in a combined optical lattice and harmonic 
trap [3|, and theoretical work investigating relaxation 
properties of the system after a quantum quench Q . 

In this Letter, we present a quantum many-body study 
of the dynamics of dark solitons formed by ultracold 
bosonic gases in one-dimensional (ID) optical lattices. 
Solitons are localized, persistent, nonlinear waves which 
appear throughout nature. Dark solitons formed by 
Bose-Einstein condensates (BECs) were first observed 
in 1999 and have lately been the topic of excit- 
ing experimental research [y, 0, Q- Owing to negli- 
gible quantum depletion out of the condensed mode, 
BECs in the mean-field limit are well-described by the 
Gross-Pitaevskii (GP) or nonlinear Schrodingcr equation 
(NLS). Such a mean-field description has proven to be 
an excellent model for describi ng g round state proper- 
ties and soliton dynamics 0, llOl . U | for experiments 
in harmonic traps. However, finite temperature [l^ and 
quantum fluctuations [l^l have both been shown to affect 
the stability of dark solitons in such continuous geome- 
tries. 

In the presence of an optical lattice, one can still at- 
tempt to describe BECs within a mean-field picture, ei- 
ther by directly employing the continuous NLS with an 
external lattice potential [3, [3| or by applying a lowest 
Bloch band tight-binding approximation to the conden- 
sate order parameter in the continuous GP equation [l6| . 



The latter procedure results in the discrete nonlinear 
Schrodinger equation (DNLS). On the other hand, the 
Bose-Hubbard Hamiltonian (BHH) is a discretiza- 
tion of the full many-body Hamiltonian that can be al- 
most perfectly realized with a system of ultracold bosons 
in an optical lattice [3. In fact, the DNLS is most 
perspicuously obtained as a mean-field approximation of 
the BHH UM- In this work, we consider quantum en- 
tangled dynamical evolution of the dark soliton — a ro- 
bust, emergent property of the corresponding mean-field 
theory (DNLS) — by performing quasiexact simulations of 
the BHH using newly available quantum algorithms [i^l . 
Hence, we directly address how many-body effects such as 
quantum fluctuations and quantum entanglement affect 
the stability of dark solitons, thereby providing a quan- 
titative measure of the applicability of mean-field theory 
in describing such dynamics. 

The Bose-Hubbard Hamiltonian describing ultracold 
bosons in a ID optical lattice reads 
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where J is the hopping strength and U is the on-site 
atom-atom interaction energy. Equation ([T]) assumes box 
boundary conditions on a lattice of M sites. The bosonic 
destruction and creation operators bk and SJ. obey the 
usual bosonic commutation relations, and fik = b^bk is 
the number operator that counts the number of bosons at 
site k. The (quasi-) ID regime is reached experimentally 
by making a 3D lattice and ramping up the lattice heights 
in the transverse directions creating an array of ID tubes 
with a sinusoidal potential in the longitudinal direction. 
Typically, each tube contains 10-100 atoms and each site 
in the tube only one or a few atoms. 

To obtain the DNLS from the BHH, one can evolve the 
bosonic destruction operator bk forward in time in the 
Heisenberg picture using Eq. ([1]): ihdtbk = [bk,H]- As- 
suming zero quantum fluctuations obtainable via a tensor 
product of on-site Glauber coherent states, bk can be re- 
placed with its expectation value ipk = {bk), resulting in 
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the DNLS to describe the condensate order parameter: 



(2) 



The solution is normahzed to A^dnls = J2k=i IV'fcP- 

In this work, we consider two methods for generat- 
ing dark sohton initial conditions in the full many-body 
BHH. In method 1, we work within a truncated, non- 
number-conserving Fock space that allows each site to 
contain at most d ~ 1 bosons. Then, via constrained 
imaginary time relaxation, we calculate the fundamen- 
tal dark soliton solution {ipk} of the DNLS and carry it 
over to Fock space by using a product of truncated co- 
herent states of the form |^) ~ (S'feLi \^k), where \zk) — 
TVd exp(-|zfcp/2)X]^^=o tItI"); ^ normalization 

factor required by our truncation and ^pk = {bk) ~ Zk- In 
method 2, we generalize the methods of density and phase 
manipulation for soliton creation [2l| to the BHH by first 
using imaginary time propagation with an external Gaus- 
sian potential to dig a density notch and then appl yin g 
an instantaneous phase gradient across the notch [22j . 
For simulation of the BHH in real and imaginary time, 
we use Vidal's time-evolving block decimation (TEBD) 
algorithm retaining Schmidt basis sets of size x- This 
method is equivalent to a time-adaptive density matrix 
renormalization group routine based on matrix product 
states [2^, 24 1 . Its accuracy depends explicitly on the 



amount of spatial entanglement in the states being sim- 
ulated [Hi. 

To characterize the quantum nature of the system we 
employ six distinct measures, (i) The quantum deple- 
tion describes the occupation of non-condensate modes. 
The natural orbitals of the system are the eigenvectors 
of the one-body density matrix (SjSi), where we denote 
the kth component of the {j 



pied natural orbital as 



l)th most highly occu- 
with corresponding occupa- 



tion Nj . The condensate wave function is the eigenvector 



of the one-body density matrix whose occupation A^o 
is largest [l^ . Depletion out of the condensate mode is 

defined as D = 1 - No/Navg, where TVavg = Z]feii(^fc) 
is the total average number, (ii) The order parameter 
(bk) maps directly onto the DNLS dependent variable 
for infinite-dimensional coherent states and has a time- 
dependent norm Nh = X^feli that measures the 
system coherence. In the thermodynamic limit in ID, it 
is well known that the local order parameter (bk) van- 
ishes; however, Bosc condensation is possible in finite 
quasi-lD systems, and products of coherent states with 
(bk) serve as a mean-field description of such con- 
densed states, (iii) The on-site expected particle num- 
ber {-hk) is the local density actually measured in exper- 
iment, (iv) The local von Neumann entropy, or entropy 
of entanglement, is an entanglement measure defined as 
SvN.k = -Tr[pk\og^{pk)] G [0,1] with pk = Trj^kP- It 
measures the entanglement of the fcth localized mode 



with the rest of the lattice 
purity 



(v) The average local im- 
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describes how far the full system is from a local pure 
state at each site, (vi) Finally, the pair correlation func- 



(J)\b^,bjbi) measures the joint probability of 



tion 5^1^ 

measuring two particles at sites i and j. 

Figure[T]depicts a characteristic simulation of standing 
quantum soliton propagation for initial conditions ob- 
tained via method 1 as described above. The parame- 
ters are vU / J = 0.35 at filling ly = A'dnls/-^ = 1 for 
M = 31 lattice sites with x = 50 and d = 7. For all 
results, we have checked for sufficient convergence in all 
numerical parameters, i.e., x, d, as well as St, the Trot- 
ter time step size used in the TEBD We refer to 
the parameter vU / J as the effective interaction strength 
because it accounts for both the atom density v and the 
interaction energy U , all scaled to J . At the initial time 
when the system is in a product of coherent states resem- 
bling a standing soliton, the depletion as calculated by 
diagonalizing the one-body density matrix is negligible 
{D K, 0.1%). The condensate wave function according to 
the Penrose-Onsager definition closely resembles the sta- 
tionary DNLS soliton solution, and this mode is occupied 
by all but -DA^avg bosons. However, unitary evolution ac- 
cording to the BHH causes an increase in occupation of 
non-solitonic orbitals giving the soliton a finite lifetime. 




FIG. 1: (color online) Three ways to describe a quantum soli- 
ton: (a) particle number density, (b) Penrose-Onsager con- 
densate wave function density, and (c) order parameter den- 
sity versus position and time during quantum evolution of a 
dark soliton initial state obtained via method 1. Shown in (d) 
is the time dependence of the average local impurity, quantum 
depletion, and order parameter norm. 
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The second most highly occupied natural orbital <b^ is 
a density maximum that fills in the soliton notch |22| . 

These quantum dynamics can be explained by noting 
that the standing dark soliton is an antisymmetric wave 
function, and there are allowed two-body scattering pro- 
cesses which preserve the symmetry of the many-body 
wave function that can deplete particles out of the soli- 
ton mode. For instance, two particles initially in the dark 
soliton orbital can scatter into a symmetric orbital which 
need not necessarily vanish at the lattice center. We see 
in Fig. [T]that such processes are energetically preferable 
over time as particles deplete into the mode (j)^^^ as well 
as orbitals of higher order. That is, our theory also treats 
higher order processes such as two atoms previously de- 
picted into a symmetric mode scattering into an antisym- 
metric mode of higher order than the soliton mode [22] . 
At time tJ/H k. 33, (j)^^'^ gains higher occupation than 
the soliton mode initially denoted This crossing in 

occupation numbers of natural orbitals is indicated by 
the discontinuity and corresponding dash-dotted white 
line in Fig. [ijb). The decay of the order parameter den- 
sity shown in Fig. [UJc) closely follows the collapse of the 
soliton structure depicted in Fig. [Ija) where the black 
dashed line corresponds to the exponential decay time rt 
of the order parameter norm Ni,. In contrast, the evolu- 
tion of the condensate order parameter according to the 
DNLS (not shown) reveals a stably propagating density 
notch, and the initial DNLS soliton solution is dynami- 
cally stable [l^]. AH in all, the behavior described above 
is general, and the qualitative picture remains when a 
harmonic trap potential is added. 

In Fig. [21 we depict the decay times tq of the soliton 
contrast for a range of vU / J values at four filling fac- 
tors: V = 0.5, 1.0, 1.5, and 2.0. The soliton contrast at a 
given time is defined as C ee {{fi^^^) - (nniid))/((nmax) + 
(niiiid))i where (rimax) is the maximum of the average 



density over all sites k and 



i) is the density at the 




FIG. 2: (color online) Quantum soliton lifetimes. The de- 
cay time Tc of the soliton contrast versus interaction strength 
vU I J at four separate filling factors v. Here, we use data cor- 
responding to dark solitons created via method 2 with M = 30 
lattice sites [S^l . All values of vU / J considered reside within 
the superfluid region of the BHH ground state phase diagram. 



lattice center. The time scale tc is the time at which C 
decays to a value of 1/2. In order to achieve higher fill- 
ing factors, we use method 2 of soliton generation with a 
number-conserving TEBD routine. Convergence of these 
results requires taking up to x = 120 and c? = 9. It is 
clear that the soliton lifetime decreases with increased 
interaction strength and increases with increased filling 
factor. 

Our approach allows a full quantum many-body char- 
acterization of the dynamics. In Fig. [ijd), wc sec that 
even though the average local impurity and quantum de- 
pletion vanish initially due to the nature of the initial 
state, these measures grow over time. This is also true 
for the local von Neumann entropy as shown in Fig. [3Kb) . 
Finally, in Fig. [3ja) we calculate g^^^ = g^j^ j,, the pair 
correlation function of site k with the center lattice site, 
i.e., the initial position of the soliton. This measure does 
not remain zero over time at all sites and is largest in 
the decaying soliton's tails; thus, the soliton decay phe- 
nomenon observed above is not the result of averaging 
a wandering, yet nondecaying, dark soliton over many 
measurements. 

In Fig. [4| we illustrate soliton-soliton collisions as de- 
scribed by both the DNLS and the BHH. The initial con- 
ditions are obtained by applying the methods of density 
and phase engineering to the DNLS using Gaussian po- 
tentials for the density engineering and hyperbolic tan- 
gent phase profiles for the phase engineering and then 
carrying the state to Fock space via truncated coherent 
states, as in method 1 for standing solitons. Figure [41(a) 
depicts the DNLS dynamics of a near-elastic soliton col- 
lision, and Fig. [HJb) shows the density during the corre- 
sponding quantum evolution. In this case, the collision 
occurs before the decoherence time Tf, and the average 
particle number closely follows the mean-field order pa- 
rameter density: the collision is still very elastic in accor- 
dance with mean-field theory. In the DNLS, the filling 
v only changes the norm of the solution for fixed vU / J] 
however, in the quantum picture, we can vary the time 
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FIG. 3: (color online) Many-body characterization of a quan- 
tum soliton: pair correlations and entanglement entropy. The 
(a) pair correlation function and (b) Local von Neumann en- 
tropy for the same simulation presented in Fig. [T] 
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FIG. 4: (color online) Quantum-induced inelasticity. Dark 
soliton collision for (a) the DNLS and the corresponding quan- 
tum evolution at filling factors (b) u = 1, (c) = 0.5, and 
(d) 1/ = 0.1 at effective interaction strength vU/J = 0.35 for 
M = 31 sites. For (b), (c), and (d), the decoherence time 
Tb [black dashed line, cf. Fig. [IJa)] is tuned to occur after, 
during, and before the collision time, respectively. 



of decoherence by varying the filling as evident in Fig. [21 
Wc show in Fig. [l{c)-(d) that the elasticity of a collision 
of two quantum solitons decreases when the decoherence 
time becomes comparable to the collision time. That is, 
the solitons interact or "stick together" for a longer time. 
In these cases, the second most dominant natural orbital 
t/)'^^) is a standing soliton, so that increased occupation 
of this mode has the effect of increasing the time over 
which the solitons collide. Subsequent filling in of the 
notch after collision can be explained by bosons occupy- 
ing even higher order modes, i.e., 4'^'^\ etc. The fact 
that the solitons can be made to stick together is another 
compelling piece of evidence that we are not observing 
diffusing solitons averaged over many measurements. 

The increase in occupation of higher order natural or- 
bitals, i.e., depletion out of the condensate wave func- 
tion, is a completely many-body effect that cannot be 
accounted for using an unperturbed mean-field theory 
such as the DNLS. However, the lowest-energy mode of 
the Bogoliubov spectrum for a stationary dark soliton 
solution of the DNLS can be an anomalous mode with a 
density maximum in the center of the lattice [isj . This 
Bogoliubov mode resembles, but is not equivalent to, the 
second order mode 0^^^ observed for the standing soli- 
ton case [iij. By simulating the BHH, a true quantum 
field theory, we are able to calculate explicitly the time 
dependence of the distribution of the natural orbitals as 



well as their exact spatial form. 

By calculating the stability times of dark solitons, 
which are stable structures in the DNLS but nonequi- 
librium states of the BHH, we are evaluating the valid- 
ity of using the DNLS to describe the system dynam- 
ics. It is thought that in the superfluid regime of the 
Bosc-Hubbard phase diagram mean-field theory should 
be applicable. However, we show that there is always a 
time at which quantum fluctuations cause such a model 
to break down, providing specific predictions in Fig. ^ 
As one tunes toward the Mott border, the stability times 
decrease, indicating an even stronger presence of quan- 
tum fluctuations and a further breakdown of mean-field 
theory. As U/J is increased past the region studied in 
Fig. [21 one eventually obtains a soliton that extends only 
over a single site, necessitating a multi-band BHH. 

In conclusion, we have constructed quantum many- 
body analogs of dark solitons and analyzed their dynam- 
ics. This was achieved by using Vidal's TEBD algorithm 
to simulate real time evolution in the Bose-Hubbard 
Hamiltonian of initial dark solitonlikc many-body states. 
Wc showed that quantum effects cause a finite soliton 
lifetime and induce an inelasticity in soliton-soliton col- 
lisions. 
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